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Abstract 

Simple analytical formulae, directly relating the experimental geome- 
try and sample orientation to the measured R(M)XS scattered intensity 
are very useful to design experiments and analyse data. Such formulae 
can be obtained by the contraction of an expression containing the polari- 
sations and crystal field tensors, and where the magnetisation vector acts 
as a rotation derivative [3]. The result of a contraction contains a scalar 
product of (rotated) polarisation vectors and the crystal field axis. As 
an example, the dipole-dipole RXS scattering amplitude by Jahn- Teller 
distortion is calculated as: 



C(e'(2z^ — — y^)e) = 2e' zze + 2e zze + ... 

= 4e'j.ez - 2eUi. - 2e'yey 

The contraction rules give rise to combinatorial algorithms which can 
be efficiently treated by computers. In this work we provide and discuss 
a concise Mathematica code along with a few example applications to 
non-centrosymmetric magnetic systems. 



1 Introduction 

A tensorial contraction method has been developed [3] to obtain analytical for- 
mulae for X-ray resonant magnetic scattering, working from the established 
formulae for RMXS amplitude in the spherical atom approximation [2]. Here 
this method will be implemented in a Mathematica code and used to investigate 
electric dipolc and quadrupole scattering in the general case of a non-spherical 
system. 
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To illustrate the capabilities of the code, we consider the case of spiral anti- 
ferromagnetic holmium. This has an HCP structure and the atoms are embed- 
ded in a local D^h symmetry environment. The magnetic field direction of the 
holmium varies helically according to the crystal layer. 

The Mathematica code executes contractions of the tensors of an incoming 
and outgoing photon with the field tensor of the crystal system, to obtain expres- 
sions for the dipolc-dipole, quadrupole-quadrupole and mixed dipole-quadrupole 
corrections to the scattering amplitude for the spherical case. 



2 Theory 

2.1 Spherical case 

We consider the interaction between matter and a photon described by a wave 
vector k and a polarisation vector e, discarding both clastic Thompson scattering 
and the spin-magnetic field interaction. Our atom is spherical [3] but perturbed 
by a magnetic exchange field. We take the angular momentum quantisation axis 
^ along the magnetic field. Since dipolar and quadrupolar terms will not mix 
due to parity conservation, the scattering amplitude is given by [3]: 

q=l q=2 

= E P^A*^\ + E ^2,,(fce);(fce), (1) 

g=-l q=-2 

where (fce)g denotes the rank 2 spherical tensor components oi k (E) e, and 
represents e'in rank 1 spherical tensor form. The scattering coefficients Fi^q 
and _F2,g can be derived from the analysis in [3]. 



2.2 Spherical case in Cartesian representation 

It is, however, more convenient to keep the vectors k and e'in Cartesian form and 
re-express Equation [T] in Cartesian space. Using a different set of coefficients 
pra^ pin (^-^^jjich -^^g -^iw define in Section and substituting q in Equation [1] 
by ^.L where L is the angular momentum operator, we rc-express Equation [1] 
as 

(n=2 \ / n=4 \ 

e'Y^iiiy^TF^'eUcUk' ®e')Y,{iixrF!,^{k®e) \ /2 (2) 
n=0 / V n=0 / 

where C signifies a sum over all possible contractions (the contractions are 
defined in Section 12. 4[) . This formulation can be used to calculate directly the 
dipolar and quadrupolar scattering coefficients for a spherical atom [3], which 
are in agreement with the formulae obtained previously by Hill & McMorrow 
[2]. The implementation of C as a computer-efficient algorithm is one of the 
objectives of this code. 
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2.3 Non-spherical case 



We use the extension of this formulation for the case of a non-spherical atom 
[3]. We represent the non-sphericity of the atomic environment by a crystal 
potential T{x, y, z) which is added to the atomic Hamiltonian, which is given as 
a superposition of spherical harmonics [3] : 

r(x,y,^) = ^fi,,T^(a:,2/,z) (3) 

l,q 

where T must also have the same symmetry as the crystal system [1]. 

The field tensor must now be included in Equation [51 The scattering ampli- 
tude for the non-spherical case is therefore given by the contraction sum 



ri=2 \ / n=4 



= C e' ^(iex)"F;"Te + C (fc' ® e') ^(iex)"Ff T(fc ® e) /2 

V n=0 / \ n=0 I 

(4) 

Later we will construct some code to evaluate i^£_fc^e',fc' in the general case 
of a field tensor in a;, ?/ and z, and then investigate the effect that this will have 
on the scattering peaks measured. 



2.4 Contractions 

First wc define mathematically what is meant by a contraction [4] and the 
function C, which represents the sum of all possible contractions between the 
tensors it is operating on, that result in a scalar. In Section [3] this will be 
implemented as a computer code. 



Single contractions A single contraction of two tensors involves taking the 
inner product between an index of the first tensor a\ 0,2 and one of the second 
tensor h\ (8) 62, thus reducing the total rank of the expression by two. So, for 
two tensors a\ ® 0,2 and h\ (g) 62, each of rank 2, the possible single contractions 
are 

(01.61)0262 {a,\.h2)a2h\ (02. 61)0162 {a2h2)a\h\ (5) 

i.e. there arc four possible single contractions of ai ® 02 with 61 (g) 62, each one 
a tensor of rank 2. 



Double contractions However, we are interested in contracting the two ten- 
sors to form a scalar. For the case of the two rank- 2 tensors a\ ® 0,2 and 61 (X) 62, 
we can execute the first single contraction and then contract the two indices 
that have not already been contracted. Working from Equation [5l we see that 
this results in a total of four possible (double) contractions: 

(01.60(02.62) (01.62) (02. 61) (02.60(01.62) (a2.62)(ai.6i) (6) 
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Therefore, the contraction sum of ai ® 02 and bi CE) 62, which is the sum of 
aU possible contractions, is 

C{A,jBki) = 2{ai.bi){a2.b2) + 2{d2.h){ai.b2) (7) 
Note the factors of 2 which arise from the repeated terms in Equation [S] 

Multiple contractions We would like to construct an iterative computer 
algorithm to take a set of tensors and repeatedly contract them until the result 
becomes a scalar. 

We begin by taking all possible single contractions (for the case of two tensors 
of rank TV, there are 27V possible single contractions; however we would also like 
to contract three tensors with each other). 

Then we take this result and take all possible single contractions on it again — 
this continues until either the expression becomes a scalar (in which case the 
contraction function is no longer called, and the result is output as the answer), 
or a tensor of rank 1 (which will be ignored by the program, as it cannot be 
contracted to a scalar). 

3 The program 

Here we describe the workings of the Mathematica notebook. The program 
input is presented in a different colour for easy differentiation from the rest of 
the document. Occasionally the program output has been stylised slightly for 
this report, but everything listed here was produced more or less directly by 
Mathematica. 

In this section we will define the contraction mechanisms for calculating the 
magnetic diffraction amplitude for the general case of a material with a known 
field tensor. 

In Section m we will then apply this mechanism for the specific case of spiral 
antifcrromagnetic holmium, and show how this can be used to calculate the 
expected intensities of the satellite diffraction peaks of Bragg order n±q. 

3.1 Tensor contraction 

First we set up a function to execute tensorial contractions (these are needed 
to calculate the scattering amplitudes for a non-spherical atom). 

Contracting two tensors First we define a function Con which will find the 
sum of all possible contractions for two tensors rank N. Each tensor is given in 
the form T[ai, d2.-.apf] (representing the tensor product ai ® a2 ^ ■■■ O oat of a 
set of vectors oi, 62, ... Sn)- For example, the contraction sum of two tensors 
ai (g) (22 and &i (8) &2 will be entered as Con[T[ai, 02], T[&i, 62]]; and is equivalent 
to, mathematically, 

C((ai ® a2)(6i ® ^^2)) = (al.6i)(a2.&2) + (ai.fe2)(a2.&i) (8) 
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i.e. each vector in the left hand tensor is contracted with each one on the right, 
thus generating all possible contractions. 

The contraction function Con for two tensors is given below. It loops round, 
gradually reducing the expression by eliminating pairs of vectors and calling 
itself again to contract the remaining expression. This creates a complicated 
stack of one function calling itself again many times with different arguments — 
to prevent infinite loops in the case of an error, in this definition Con will not 
call itself directly but will call a temporary (as yet undefined) function ConTmp 
which can later be set equal to Con. 

First we state that the function Con and our tensor product operator T are 
both orderless, so it does not matter in which order their arguments (vectors 
and tensors) are given. 

Attributes[Con] = {Orderless}; 
Attributes[T] = {Orderless}; 

Now we create the contraction algorithm. The function takes two tensors 
and returns the sum of all possible single contractions between them. 

Con[x__, 0]:=0; 

Con[T[a_, al ],T[bl ]]:=Module[{res}, res = 0; 

n = Length[List[bl]]; 

Do[ 

res — 

res + seal [a/, {x ^ ei, y ^ e2, z — > 63} , 
List[bl][[i]]/. {x ^ ei,y ^ e2,z eg}] 
X ConTmp[T[al],Delete[T[bl], i]] 
,{i,l,ii} 

]; 

res 



where we introduce the variables ei , 62 and 63 to represent the basis vectors in 
the X, y and z directions (this will be useful later for when we will need to take 
scalar products). 

Some other definitions are necessary to enable the program to stop looping 
when it has finished the contractions. The contraction of a single tensor must 
evaluate to zero; to contract anything with 1, we can ignore the 1; and the 
contraction of nothing must evaluate to 1 . 

Con[T[c___]]:=0; 

Con[l, t ]:=Con[t]; 

Con[]:=l; 

TD = 1; 
T[{}] = 1; 

T[{x..}]=T[x]; 
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Contracting three tensors Now a contraction function has been defined 
for two tensors, we can start to think about C for three tensors. This is more 
comphcated, as the sum of possible contractions is much greater. The function 
works in the same way as the Con function for the case of two tensors. cHmi- 
nating vector pairs one by one and caUing itself (via ConTmp) repeatedly. Once 
one of the three tensors has been fully contracted with vector elements from the 
other two, the expression will have been reduced to a sum of contractions of two 
tensors — which is a much simpler problem and has already been defined above. 

Con[T[al___], T[bl___], T[cl___]]:= 
Module[{res, m, n}, res — 0; 

m = Length[T[bl]]; 

n = Length[T[cl]]; 

Do[ 
res = 

res + seal [List[al] [[1]]/. {x ^ ei, y ^ e2, z ^ ea} , 

List[bl][[i]]/. {x ei,y ^ e2,z ^ eg}] 

X ConTmp[Delete[T[al], 1], Delete[T[bl], i], T[cl]] 
,{i,l,m} 
]; Do[ 
res = 

res + seal [List[al] [[1]]/. {x — > ei, y — > e2, z — > ea} , 
List[cl][[i]]/. {x ei,y e2,z eg}] 
X ConTmp[Delete[T[al], 1], T[bl], Delete[T[cl], i]] 
,{i,l,n} 

]; 

res 

]; 

As a temporary measure, we also make all callings of ConTmp equivalent to 
calling the function Con itself (thus creating a stack) — in the case of errors or 
infinite loops, this line can be removed. 

ConTmp — Con 

3.2 Extending the contraction functions to work with poly- 
nomials 

The field tensor in polynomial representation Now we need to make 
a function which will contract a Cartesian tensor in the form of a polynomial 
in a;, y and z. This is needed for scattering amplitude calculations, where the 
crystal field tensor is typically expressed in a polynomial form T = f{x,y,z). 
For example, the quadrupole-quadrupole scattering formula for the particular 
case of holmium has the form: 
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F^:L.'^k' = k'){t2i2z^ - x^~ y^) ± tsix"" ~ ixy^W ® k)) (9) 



Wc define a new function C{A, f(x, y, z), C) which re-expresses a contraction 
of tensors A and C and the polynomial / in terms of sums and products of 
contractions of three tensors which can be executed by the function Con. 

Reducing the polynomial via linearity We want to be able to contract 
some quite complicated polynomials, so it is useful to first reduce the polyno- 
mials to a more manageable form using some basic mathematical properties of 
the function C. For example, the contraction sum operation is linear, so 

C{Aif{x, y, z) + g{x, y, z))B) = C(A/(x, y, z)B) + C{Ag[x, y, z)B) (10) 

So we can start constructing the function C by creating a Rul^ll that any 
polynomial to be contracted must first be split into terms and the final contrac- 
tion will be a sum of contractions of these terms: 

C[{al..},Plus[x.,y_.],{cl..}]:= 
Module[{i}, xy = Expand[x + y]; 
Sum[C[{al}, xy[[i]], {cl}], {i, 1, Length[xy]}]] 

Removing prefactors from the individual terms Now we need to find an 
efficient way of taking each term and extracting the xyz dependent components 
for contraction, while leaving numbers and constants untouched by the Con 
function. Here we use another property of the linearity of C: 

C(A(fc/(.T, y, z))B) = fcC(A/(.T, y, z)B) (11) 

So first we separate out any occurrences of xyz from the terms of the poly- 
nomial (which has already been split up into a sum of parts): 

C[{al..}, 

Times[prefactor_/;FreeQ[FullForm[pref actor], x]/; 
FreeQ[FullForm[pref actor], y]/; 
FreeQ[FullForm[pref actor], z], xyz_], {cl__}]:= 
Times[pref actor]C[{al}, Timesfxyz], {cl}]; 

Converting a polynomial term to a tensor Now that individual terms in 
xyz have been separated, they can be converted into a tensorial form (such as 
T[x,x,y, z] in Mathematica) that can be used directly by the function Con. 

^In Mathematica, a Rule defining C means that every time C is referred to in subsequent 
code, the code within the Rule is invoked. Parameters on the left side of a definition of a Rule 
must be followed by underscores (e.g. A_), to show that they can take any value. 
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Expressing x°'y^z'^ in tensorial form First we make a Rule that C{A{x°'y^ z^)B) 
will be executed as Con[T[A], T[B], T[C]] where B is a sequence of the variables x, 
y and z. For example, we want the expression 

C[T[A], (x^yz),T[C]] = C{A{x'^yz)B) {mathematical notation) (12) 

to be converted to 

Con[T[A],T[x,x,y,z],T[C]] (13) 

and then contracted in the usual way by the function Con. The code to execute 
this is as follows: 

C[{al__}, Times[xyzl_, xyz2_] , {cl_}] := 
Module[{i, power, n, xyz}, 
xyz:=List[xyzl, xyz2]/.x_P°'"'''-- >Table[x, {i, 1, power}]; 
n = Length[xyz] ; 

Con[T[al],T[Flatten[Table[xyz[[i]],{i, 1, n}]]], T[cl]] 

] 

Expressing x" in tensorial form We also need a Rule to interpret powers 
of a single variable, such as C{A{x")B), since this case is not covered by the 
above Rule (although the interpretation is much the same). 

C[{al__}, Power [xyz_, power_] , {cl_}] := 
Module [{i,n}, 
Con[T[al],T[Table[xyz, {i, 1, power}]], T[cl]] 

] 

3.3 Defining the scalar product symbolically 

The contraction function defined so far has produced results in terms of the func- 
tion seal [a, 6], as yet undefined. We must now define seal so that a product of 
a Cartesian vector with one of the basis vectors x, y and z will be re-expressed 
as the correct component of that vector (at least for the moment). 

Clear[scal, ^]; 

Attributes[scal] — {Orderless}; 

seal [vector., eij :=veetor.{x, y, z}[[i]]; 

4 Application: spiral antiferromagnetic holmium 

The contraction mechanism which has been defined will work in general on any 
kind of polynomial tensor. Now we present an example of its application, for 
the specific case of spiral antiferromagnetic holmium 
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Figure 1: The geometry of the holmium system. Since r is different for different 
crystal planes, Zm has a spiral dependence on z. 

4.1 Defining the geometry of the system 

Defining frame X First, we define a simple Cartesian co-ordinate system. 
Initially we would like to work in the frame X with co-ordinates (x, y, z) oriented 
around the crystal geometry (z points upwards through the crystal planes), as 
shown in Figure [47T] 

The magnetic field direction Zm of spiral antiferromagnetic holmium varies 
helically according to the crystal layering, as 

^m(T) = (cost, sin T,0) (14) 

where the angle t = ^^j^ depends on the height z [d is the interplanar spacing 
and q is the antiferromagnetic wavevector). In Mathematica code this is ex- 
pressed as (see Figure |4?T|) : 

z„,[t.] = {Cos[r],Sin[T],0}; 

The incoming/outcoming photons Now we define the characteristics in 
frame X of an incoming photon, at angle 9 from the crystal planes. Its wavevec- 
tor is 

kx[6l.] = {Cos[6l],0,-Sin[6l]}; 

and its polarisation vectors in X are (with the indices and 1 representing 
polarisation orientations a and tt respectively): 
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exo[0_] = {0,1,0}; 
exi[6i_] = {Sin[6i],0,Cos[6']}; 
vector^ :=vectoro; 
vector_7r:=vectori; 

4.2 The magnetic co-ordinate system 

Defining M It is more convenient for our purposes to introduce a new co- 
ordinate system M with basis vectors {xm,ym, Zm) here, constructed from the 
(rotating) magnetic moment Zm- We set Xm{T) to the most convenient possi- 
biUty: 

x,[r.] = {0,0,1}; 

And we can now derive ymi^) from 2;,„(t) and z„i{t): 

ym[T_] = Xm[T] X Zm[T]; 

So our new basis is 



/ — sin T cos T \ 

{xm\ym\Zm) = COST sin T (15) 

' 1 ' 



Re-expressing wave parameters in the magnetic co-ordinate system 

Now we can redefine k and e^.tr in the frame M from tlicir definitions kx and 
ex^ o- in frame X (the subscript pol stands for either tt or cr): 

k[0_,T_]:={kx[^?].x„[T],kx[0].y„[r],kx[0].z,[r]}; 

epoi.[6'_, r_]:={eXpoi[6'].x„[r],eXpoi[6'].ym[r], eXpoi[6'].Zm[r]} ; 

So the wave vector in this new basis is 

- sin(6l) 

fc(6l, r) = ( - cos(6') sin(T) | (16) 
cos{9) cos(t) 



And the polarisation vectors are 











cos(r) 






V sin(T) j 





cos(6') 

- sin(6') sin(r) | (17) 
cos(r) sin(0) 
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The reflected wave We define the reflected wave parameters k' , by chang- 
ing the sign of 6: 



Cpoi [^-7T-]:=epoi[-6',T]; 

4.3 Applying the contraction mechanism for the Holmium 
case 

Using our expression for the Holmium field tensor (ignoring the ± alternation), 

T = t2{2z^ -x^ - y^) + t3(a;3 - Sxy^) (18) 

we can use the functions C and Con to derive some properties of the diffrac- 
tion pattern, that arise from the non-sphericity of the holmium atom. We are 
interested in particular in effects that arc unique for the non-spherical case — 
corrections to the spherical scattering formulae arc what wc arc after. We will 
analyse separately the dipole-dipole, quadrupole-quadrupole and mixed terms, 
as the whole expression is rather complicated. 

Calculating the dipole-dipole scattering correction First we look at the 
simplest case: the dipole-dipole scattering correction. Here we only consider the 
contraction of e' and e (the outgoing and incoming photon polarisation tensors 
respectively) with T, ignoring k' and k: that is, we ignore effects of the finiteness 
of the speed of light. We will attempt to find this function's proportionality (i.e. 
we are not interested in prefactors for now). 

The contraction is (taking only the first order term in the operator (i^x)"): 

C(e'T(j^ X £)) - [e' ^ e] (19) 

To encode this contraction in Mathematica, we ignore for now the symmetry 
term — [e' e], and enter: 

C [{e'}, t2 (2z2 - x2 - y2) + t3 (x^ - Sxy^) , {i^ x e}] / X2 1 

- 2e'.x{i^ X e).x - 2e'.y(i^ x e).y + 4e'.z(i^ x e).z (20) 

Condensing the spherical components into a single term gives (still ignoring 
the symmetry term) 

Simplif y[%, (e'.x(iC x e).x + e'.y(i^ x e).y + e'.z(i^ x e).z) 
e'.(i^ X e)] 

-2e'.(z^ x e)-h6e'.z(i^ X e).z (21) 
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Subtracting the symmetry term [e' <-!■ e], and removing constant prefactors, 
we see that the dipole-dipole scattering correction is proportional to: 



Simplif y[%/2 — syinm]//TraditionalForm 



symm — e'.(i^ x e) + 3e'.z(i^ x e).z 



(22) 



This is Equation 29 in [3]. The e'.(i^ x e) part merges with the term that 
appears in the spherical formulae, but the other term adds complexity to the 
amplitude dependence on experimental geometry. 

Calculating the quadrupole-quadrupole scattering correction Math- 
ematically, the contraction we are interested in is 



C((fc' ® e')T{i^x){k (g) e)) - [fc', e' <^ k, e] 
= [C((fc' ® e')T{{iC X fc) ® e)) + C((fc' (g, e')T{k ® {iS. x e)))] (23) 



In Mathematica, this is encoded as (ignoring the {i^ x fc, e} term, and also 
the symmetry terms, for now): 

C [{k', e'}, t2 (2z2 - x2 - y2) + t3 (x^ - 3xy2) , {k, x e}] /. 
t2 ^ 1 

which gives output 

-2e'.x(i^ X e).x scal[A:, k'] - 2e'.y{i^ x e).y scal[fc, k']- 
2fc'.x(i^ X e).x scal[A;,e'] -2k'.y{i^ x e).y scal[A:,e']- 

2fc.x(e'.x scal[fc',i^ x e] + fc'.x scal[e',i^ x e])— , , 

2fc.y(e'.y scal[A:', x e] + k' .y scal[e', il x e])+ 
2(2e'.z(i^ X e).z scal[A:,fc'] + 2fc'.z(i^ x e).z scal[fc, e'] + 
2fc.z(e'.z scal[fc', x e] + fc'.z scal[e', x e])) 

Condensing the spherical components into a single term, 
Simplify[%, 

{e'.x(i^ X e).x + e'.y(iC x e).y + e'.z(i^ x e).z == e'.(i^ x e), 
k'.x(i^ x e).x + k'.y(i^ x e).y + k'.z(i^ x e).z == k'.(i^ x e), 
k.xe'.x + k.ye'.y + k.ze'.z == k.e', 
k.Xk'.x + k.yk'.y + k.zk'.z == k.k'}] 



-2{e'.{i£, X e) scal[A:, k'] + fc'.(i^ x e) scal[A:, e']- 
3fc'.z(i^ X e).z scal[A:, e'] + k.e' scal[fc', x e] — 
3e'.z((i^ x e).z scal[fc. A:'] + k.z scal[A:',i^ x e]) + 
k.k' scal[e',i^ x e] —Sk.zk'.z scal[e',i^ x e]) 



(25) 



and now taking only the non-spherical terms in 
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Simplify[%/.{e'.(i^ x e) 0,k'.(i^ x e) ^ 0,k.e' ^ 0,k.k' 0}] 



6(e'.z((i^ X e).z scal[fc, /c'] + k.z scal[A:',i^ x e]) + 
fc'.z((z^ X e).z scal[fc, e'] + fc.z scal[e', x e])) 



(26) 



Adding the {i^ x fc,e} term back in now (by exchanging k and x e and 
adding the result on to the original formula), and rewriting seal in the dot- 
product notation, we find that the quadrupole-quadrupole scattering correction 
is proportional to: 



which is equivalent to Equation 30 in [3] . 

Calculating the mixed (quadrupole-dipole) scattering correction For 

the quadrupole-quadrupole and dipole-dipole scattering corrections calculated 
above, only the t2 (quadratic order) terms of T had any effect. The alternating 
term ±t3{x^ — 3xy^) contributes to the quadrupole-dipole and dipolc-quadrupole 
scattering factor corrections, due to its cubic order. 

The mixed scattering correction is (ignoring for now the extra terms in 
[e <-> k] and [k fc', e <-> e']): 

Expand [C [{e'}, t2 {2z^ - x^ - y^) + t3 (x^ - Sxy^) , {i^ x k, e}] /. 
{t2 ^ 0,t3 — > 1}] + syiiim//TraditionalForm 



where symm represents [e <-> k'] + [k ^ k' ,e ^ e']]. 

After rearrangement according to the rules of the scalar triple product, this 
is equivalent to Equation 31 in [3]. This is the scattering factor contribution 
from the alternating term ±^3(0;'^ — ixy^), which is a dipole-quadrupole term. 

4.4 The physical meaning of the contraction results 

Extra peaks arising from asphericity of T The asphcricity of the holmium 
atom results in extra peaks of Bragg order 2n ± g (where n g N and q is the 
antiferromagnetic wavevector) appearing alongside the peaks of order 2n that 
one would normally expect to observe. Here we will examine the magnetic 
contribution to the scattering amplitude of the peak appearing at 2n + q, for 
four polarisation channels (ct ^ ct, tt ^ ct, tr — > tt and tt ^ tt). 



Simplif y[% + (%/.k tmp/.e -> k/.tmp -> ()]/&/. 
scal[vecl_, vec2_] vecl.vec2//TraditionalForm 



e'.z{k.zk'.{i$, x e)+ 
k'.{ii X k)t.z + k'.e{i(, x fc).z + k.k'{iS, x t).z)+ 
fc'.z(e.ze'.(i^ x fc) + k.ze'.{i^ x e) + 
e.e'(i^ X k).z + k.e'{iS, x e).z) 



(27) 



symm + 6e.xe'.x(i^ x A:).x — 6e.ye'.y(i^ x fc).x 
6e.ye'.x(i^ x k).y — 6e.xe'.y(z^ x k).y 



(28) 
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From [3], the contribution to the scattering amplitude at this peak, for a 
given field tensor T, is given by 

=c(e'T^(*^x)"^^{"ej+c((fc'0e')r ^ {i^xT F!^\k ® M / 2 

\ n = l / \ n=l,3 / 

(29) 

which originates from Equation^ where only the terms in odd n are considered. 
The F prefactors are defined as follows: 



Fl'o — Fio; 

/ _ (Fii-Fi-i) . 

i' 1 1 — 2 ' 

Tr^l (2Fio— Fii— Fi_i) . 

1 2 — 2 ' 

F2'o — F20; 

(F2-2^F22 + 8F2i— 8F2-1) . 

F2 1 - j2 ' 

n/ (I6F21 + I6F2-1— F2-2— F22— 3OF20) . 

F2 2 - 24 ' 

n/ (F22— F2-2+2F2-1 — 2F21) . 

3 - 12 ' 

■no/ (6F20— F22 — 2F2-2— 4F21— 4F2-1) . 

^2 4 - 24 ' 



where the Fi^q and F2,g originate from the treatment of the Holmium atom 
as the spherical case with a magnetic field perturbation, although we are not 
concerned with their definition here. 



Redefining the scalar product We now want to create another definition of 
seal, so that it decomposes scal[a, 6] into axhx + ayby + azhz — this is necessary 
as we will be using the definitions from Equations [16] and [17] to calculate the r 
and 9 dependence of the scattering amplitude. 

First we define the scalar product function for various combinations of vec- 
tors, basis vectors and cross products with ^. We are making a new definition 
of seal, and so we remove the previous one. 

Clear[scal, £]; 

Attributes[scal] = {Qrderless}; 

Now we define the scalar product of various expressions involving a cross- 
product with the magnetic basis vector ^, and any of the basis vectors ei, 6*2 
or 63. In our co-ordinate system, ^ = (0,0,1) = £?3, so seal has been defined 
accordingly. 

seal [i^ X (i^ x (i^ x vector.)), eij := 
i {— vectory, vectorx, 0} [[i]]; 

seal [i^ X (i^ X veetor_), ei_] :~ jveetorx, vectory. O} [[i]]; 
seal [i^ X veetor_, e^J :=i {— veetory, veetorx, O} [[i]]; 
seal [veetor_, eiJ :=veetor{x,y,z}[[i]]; 
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Now we make a similar definition for scalar products with other vectors, such 
as k or e. 

scal[i^ X (i^ X X vecl_)), vec2_]:= 
i (— veclyVec2x + veclxVec2y) ; 

scal[i^ X (i^ X vecl_), vec2_]:= (veclxVec2x + veclyVec2y) ; 
scal[i^ X vecl_, vec2_]:=i (— veclyVec2x + veclxVec2y) ; 
scal[vecl_, vec2_]:= 

Sum [vecl{x,y,z}[[i]]vec2{x,y,z}[[i]], {i, 1, 3}] ; 

Considering the expected form of the solution Knowing that for the 
case of this particular peak the magnetic contribution to the scattering ampli- 
tude will have terms in i(Fi,i — -Fi^-i), 1(^2.1 — F2.-1) and i{F2.2 ~ ^2.-2), we 
can initialise a small array, magterm, to store these three terms. 

magterm — {0, 0, 0}; 

Calculating the first part of the scattering amplitude expression Now 

we will try and find the first term, which will have a prefactor oii{Fi^i —Fi^-i)). 
This term derives from the dipole-dipole contraction: 

magterm[[l]] ~ 
Simplif y[ 

Fl'iC [{e'}, t2 (2z2 - x2 - y2) + t3 (x^ 
t2 ^ l/.Fn -Fi_i ^ \] 

After simplification, the result is 

<^v^'x - ^X^'y (30) 
i.e. there is a contribution to the scattering amplitude of 

I (Fi.i - Fi,_i) {eye', - e.^^) (31) 

This will be evaluated later on in more detail from the definitions of e, e', k 
and k' in Equations [16] and 1171 ^-iid re-expressed in terms of the incident angle 
9 and the antiferromagnetic parameters of the holmium system. 

Calculating the second and third parts of the scattering amplitude 
expression We can now calculate the terms in i{F2.i — F2.-1) and 1(^2,2 — 
^2,-2)- These originate from the contraction 

(f^^C {{k' ® e')T{iix)(k' <g) e')) + C ((P ® ^)T{iix)^k' <E) e'))) /2 (32) 



-3xy2),{iex6}]/. 
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where an operation of i^x on a tensor is expanded binomially as a sum of 
operations on the vector components of the tensor. In Mathcmatica, we calcu- 
late 

basef ormagterms2and3 = 
FullSimplify[ 
F2'i (C [{k', e'}, t2 (2z2 - - y^) + t3 (x^ - Sxy^) , 
{k,iexe}] + 

C [{k', e'}, t2 (2z2 - x2 - y2) + t3 (x^ - Sxy^) , 

{iexk,e}])+ 

F2'3 

(C [{k', e'}, t2 (2z2 - x2 - y2) + t3 (x^ - Sxy^) , 
{k, if x(iex(iCx £))}] + 

3C [{k', e'}, t2 (2z2 - x2 - y^) + t3 (x^ - 3xy2) , 
{iCxk,i^ X (iexe)}] + 

3C [{k', e'}, t2 (2z2 - x2 - y^) + t3 (x^ - 3xy2) , 
{iCx (iexk),iexe}] + 

C [{k', e'}, t2 (2z2 - x2 - y2) + t3 (x^ - 3xy2) , 
{iCx (i^x (iexk)),e}])/.t2^1]/2 

which gives the rather complicated result 



(fc. (fc'. {eye', - e^e'y) + {-k'ye, + k'^ey) (i^2.-i - ^2,1) + 

kx {k',e,e'y{-F2,-i+F2,i) + 
k'y (e,e; (-F2,_i + F2,i) + 2 + e,,e;) (i^2,-2 - ^2,2)) + 

2fc', {eye', ~ e.,e'y) {-¥2,-2 + F2.2)) + (33) 
A:j,(fc',e,e^ (^2.-1-^2,1) + 
k'x (e.e; (F2,„i - F2,i) - 2 (e,e^ + e^^e^) (^2,-2 - ^^2.2)) + 
2k' y {eye', - e.,e'y) {-F2.-2 + F2.2))) 

which can be used to calculate the second and third terms in F . 



Calculating the second part of the scattering ampHtude expression 

First we use Equation 1331 to calculate the scattering amplitude term in i(i^2,i — 
^2,-1)- This involves taking only the coefficients of 1^2,1 (by symmetry, the 
— ii^2,-i terms are identical), and discarding the other terms. 

magterm[[2]] = 
FullSimplify[ 

basef oriiiagterms2and3/.F2i ■J/.F2-1 O/.F22 O/.F2-2 ^ O] 



\ {k, {k', {-eye', + e,e'y) + {k'ye., - k',ey) + 
{-ky {k'^e', + k',e'J + k, {k' ,e'y + fc'.e'J)) 
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Calculating the third part of the scattering amplitude expression 

And now wc find the term in i(i^2,2 — ^2,-2), by the same method. 



magteriii[[3]] — 
FullSimplify[ 

basef oriiiagterms2and3/.F2i O/.F2-1 O/.F22 y/-F2-2 ^ O] 

ky (k'y (e^e^ — ^^^y) + k'x {^x^'x + ^y^y)) ~ /oc\ 
kx {k'x + e.e;) + k'y (e.e^ + e^e^)) ^""^^ 



Reducing the scattering amplitude expression The three terms in F 
that we have just calculated are: 

magterm//TableForm 



ey4 - <^xe'y (36) 
i [k. {k', {-Eye'^ + e^e'y) + (fc'.e, - fc'.e^) e'^) 

(fc'y (fyei: — ^x^'y) + fc'a; (e^e^c + ^y^'y)) /-qo-i 

-fc, (fc', (-e.ye;, + e,e;)+fc'j,(e,ei + e,ye;)) ^ ^ 

We introduce a shorthand notation and re-express the tensor products. 

i {Fii — Fi_i, F21 — F2-1, F22 — F2-2} 
X FullSimplif y[magterm, 

|6xky + 6ykx = = kexy; f^X^z ^" ^z^x^^kexz, Cykz + ezky^^kCyz, 

^x-'^ y ^ ^y-'^ X — =k e xy^ ^x"'^ z ^" ^2)^ X — — k 6 xzi 

^yk z ^ f^z-^ y=^=^k e yz, k x^x k yey=~k e x2my2 7 
kxEx ^ kyey==kex2my2}] //TableForiii//TraditionalForm 

So the reduced (simplified) form of the terms is 



*(ey4-e.4)(Fi,i-J^i,-i) (39) 
(F2,i - F2^_i) (fce,,,fc'e'j,,, - /cej,,,fc'e',,,) (40) 

« (^2.2 - F2-2) {k'e'.^2my2kex,y - ke^2my2k'e'x,y) (41) 

Finding the scattering amplitude matrices Now we will evaluate the 
expression for the scattering amplitude. Since e and e' can both be in either 
the TT or the a polarisation, we construct a 2 x 2 matrix to represent the four 
combinations of e ^ e' in terms of 6 and r. 
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First wc evaluate the matrix F ~ { 1 or rather its three 

constituent terms which originate from the scattering amplitude expression cal- 
culated above. 



Fmatrix:— 
Table [ 
Simplif y[ 

magterm[[i]]/. {vector_xyz_ 

{7ra',0, l},{7ra,0, 1}]; 
The three terms are 



vector[6',T][[xyz]]} /. 
l,y^2,z^3}],{i,l,3}, 



Table[Fmatrix[[i]]//MatrixForm, {i, 1, 3}]//TableForm// 
TraditionalForm 



-c(6')c(t) 
c(6i)c(r) -2c(e')s(6')s(T) 

/ -ic(2r)s(2^?)s(r) ^ 

ic(6')c(T)(-2c(2e') + 3c{2{0 - r)) - 2c(2r) + 3c(2(6l + t)) 
ic(6')c(T)(2c(26') - 3c(2(6i - t)) + 2c(2t) - 3c(2(6l + r)) - 2) 



v 



ic2(r)s(40)s(r) 



(42) 



/ -2c(^)c2(t)s(0)s(t) 

{0)c{t) {sH9) (2s2(r) 



\ 



c(0)c(r)((c(2r)-2)s2(0)+c2(^?)s2(r)) 
V |s(4(?)(s(3r)-7s(r)) J 

where the last two 4x4 matrices have been condensed into a columnar form, 
and the abbreviations s and c stand for sin and cos respectively. 



Fourier analysis of the F matrix Since we are interested in the +q satel- 
lite, we need to find the e*"^ Fourier component of the matrix F. The following 
code breaks F into e™"^ components and extracts the first Fourier coefficients. 

Fourierl = 
Table[ 
TrigReduce[ 
PadRight [Coef f icientList[ 
Expand[TrigExpand[Fmatrix[[i, varl, var2]]]/. 
Sin[n..T] ^ ^ (e^"^ - e-^"^) /. 

Cos[n r] ^ I (e^"^ + e-^"^)] /.e^^ ^ TERM, TERM] ,2] [[ 
2]]],{i,l,3},{varl,l,2},{var2,l,2}]; 

So the coefficient of the scattering amplitude matrix multiplying the Fourier 
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component is: 



Do [Print [" +",i{Fn -Fi_i,F2i -F2-i,F22-F2-2}[[i]]// 
TraditionalForm, Fourierl[[i]]/ /MatrixForm// 
TraditionalForm], {i, 1, 3}] 

+^(^^'^-^^-^)(, ^(7c(3^)-3c(^)) ^zs(4e) J 

The above expression is the scattering amplitude matrix for the n + q satel- 
hte diffraction peak of holmiun:|^. This can be used to obtain theoretical predic- 
tions for the amplitudes of the various scattering peaks for photons of different 
energies — which can then be compared with measurements made at the ESRF. 



5 Conclusions 

We have taken the contraction method established in [3] and implemented it 
as an efficient computer algorithm in a Mathematica program. The program 
was designed to take a general (polynomial) field tensor and contract it with 
various combinations of k and e with the magnetic vector ^ acting as a rotation 
derivative. The program output scattering amplitude expressions in terms of 
components of k and e. 

An example of the program operation was presented here for the case of 
spiral antiferromagnetic holmium, where the vector ^ rotates round through 
the crystal layering structure with an antiferromagnetic wavevector q. Using 
the non-spherical field tensor for this material, the dipole-dipole, quadrupole- 
quadrupole and mixed dipole-quadrupole corrections to the scattering factors 
which were obtained for the case of a spherical atom [2] were calculated. 

The contraction method was then used to find the amplitude of the diffrac- 
tion peak with Bragg number 2n ± q, for the a' ^ cr, c' ^ tt, tt' ^ cr and 
tt' ^ TT polarisations, via a Fourier analysis of the scattering amplitude func- 
tion. These scattering amplitudes could then be used to predict the intensities 
of the satellite peaks observed in RMXS experiments on Holmium. 

Further work This program can be developed further to allow for more com- 
plex forms of field tensors, such as those involving cross products with the 

^This can be shown by considering components of a wave incident at angle 9 being reflected 
off adjacent crystal layers zi, Z2 with scattering coefficients e"^ and e''^^. The electric fields 
of the two reflected waves will be given by e*'^ie*'°^ and ^^^2 g^^k(x+2dsln9) ^ The condition for 
these to be in phase leads to 2d sin = qX, a modified Bragg diffraction condition. 
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photon tensors (see Appendix |^ . The rehabihty of the contraction method [3] 
in predicting scattering peaks for other materials is yet to be tested. 

The extension of the program to cope with contractions of more than three 
tensors is relatively trivial, involving only the addition of a few more rules to 
reduce the number of tensors in the same way that the present program reduces 
it from three to two. 

A more complicated problem would be to contract a set of tensors with 
restrictions on the possible contractions, such as in Equation 25 of [3]. 

6 Acknowledgements 

I am grateful to Alcssandro Mirone for his advice on the programming and for 
his patience in explaining the theoretical details of his paper [3] . 

A Extension of the program 
A.l Tensor contraction 

Two tensors First we define the contraction sum for two tensors of any rank. 

Attributes[Con] — {Orderless}; 
Attributes[T] — {Orderless}; 
Con[x__, 0]:=0; 

Con[T[a_, al ],T[bl ]]:=Module[{res}, res = 0; 

n = Length[List[bl]]; 
Do[ 

res = res + seal [a/, {x ^ ei, y ^ e2, z ^ es} , List[bl] [[i]]/. {x ei, y e2, 
ConTmp[T[al],Delete[T[bl], i]] 

]; 

res 

]; 

Con[T[c___]]:=0; 

Con[l, t ]:=Con[t]; 

Con[]:=l; 

TD = 1; 
T[{}] = 1; 
T[{x..}]=T[x]; 



Three tensors Now we define the contraction sum for three tensors (by first 
eliminating one of the three, and then reverting to the two-tensor contraction). 
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Coii[T[al ],T[bl ],T[cl ]]:=Module[{res, m, n}, res = 0; 

m = Length[T[bl]]; 
n = Length[T[cl]]; 
Do[ 
res = 

res + seal [List[al] [[1]]/. {x ^ ei, y ^ e2, z ^ ea} , 
List[bl][[i]]/. {x ^ ei,y -i- e2,z — i- eg}] 
ConTmp[Delete[T[al], 1], Delete[T[bl], i],T[cl]] 
,{i, l,m} 

]; 

Do[ 
res = 

res + seal [List[al] [[1]]/. {x ^ ei, y ^ e2, z — !■ es} , 
List[cl][[i]]/. {x ^ ei,y ^ e2,z ea}] 
ConTmp[Delete[T[al], 1], T[bl], Delete[T[cl], i]] 

]; 

res 

]; 

As a temporary measure (for debugging), wc make all ConTmp be executed 
as Con. 

ConTmp — Con; 

A. 2 Contracting polynomials with cross products 

We will construct several functions, all defined by rules, along which the contrac- 
tion is passed. Each function executes its own simplifications and contractions, 
and then invokes the next one in the sequence. 

Expanding brackets The function C rewrites the polynomial in expanded 
form, and passes the contraction to C2, the second contracting function. 

C[al_, anything., cl_]:=C2[al, Expaiid[anything] , cl]; 

Expanding contractions of sums If we are a contracting a sum of terms, 
then C3 splits it into a sum of contractions... 

C2[al_, Plus[tliingl_, morethings_], cl_]:= 

Sum[C2[al, Join[List[thingl], List[iiiorethings]][[i]], cl], 

{i, 1, Length[Join[List[tliingl], List[morethings]]]}]; 

...until we we have nothing left to expand into sums: then it passes the con- 
traction on to C3, the third contraction function 
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C2[al_, Sum[onething_], cl_]:=C3[al, onething, cl]; 

Taking sums outside x operator Now C3 expands the sums within cross 

products into sums of contractions 

C3[al_, 

Times[pref actors , vecl_ x Plus[vec2a_, vec2b__]], cl_]:= 

Times@@Complement [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] 
Times [ 

Sum [C3 [al, Times@@Intersection [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] 

vecl X (vec2a + vec2b)[[i]],cl],{i, l,Length[vec2a + vec2b]}]]; 

C3[al_, 

Times[pref actors , Plus[vecla_, veclb__] x vec2_], cl_]:= 

Times@@Complement [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] 
Times[ 

Sum [C3 [al, Times@@Intersection [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] 
(vecla + veclb)[[i]] x vec2, cl], {i, 1, Length[vecla + veclb]}]]; 

Taking powers outside x operator Now C3 can also expand powers within 
cross products, and re- write the expression as a sum of cross products with each 
of the variable on the right of the cross-product operator. 

C3[al_, Times [pref actors , Cross [a_. Power [xyzl_, power.]]], cl_]:= 

Module[{i, xyz}, xyz = Table[xyzl, {i, 1, power}]; 
Times(Q)@Complement [{pref actors}, {x, y, z, x, y, z, x", y-, z"}] 
Sum[ 
C3[al, 

Times [Times@@Intersection[{pref actors}, 

{x, y, z, x, y, z, x-, y-, z-}] , Times@@Delete[xyz, i], cross[a, xyz[[i]]]] , cl] , 
{i, 1, Length[xyz]}] 

] 

Taking products outside x operator C3 also expands products within 
cross products, similarly re-expressing the polynomial as a sum of cross product 
operations on each of the variables on the right of the x sign. 

C3[al_, Times [pref actors , Cross [a_. Times [xyz 1_, xyz2__]]], cl_]:= 

Module[{i, xyz}, 

xyz — Flatten[Join[List[xyzl], List[xyz2]]/. 
{anything_P°"^^- Table [anything, {i, 1, power}]}] ; 
Times@@Complement [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] 
Sum[ 
C3[al, 

Times [Times@@Inter sect ion[{pref actors}, 

{x, y, z, x, y, z, x", y-, z-}] , Times@@Delete[xyz, i], cross[a, xyz[[i]]]] , cl] , 
{i, l,Length[xyz]}] 
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] 



Introducing the cross operator We now re-express any simple cross prod- 
uct witliin a contraction as a function of the cross product operator. 

C3[al_, Times [pref actors , Cross[a_, b_]], cl_]:= 

Times@@Complement [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] 

C3 [al, Times@@Intersection [{pref actors}, {x, y, z, x, y, z, x", y-, z-}] cross[a, b], 



Cross[(a/.{x ^ x,y ^ y,z ^ z}/-{x ^ {1,0,0}, y ^ {0, l,0},z ^ {0,0, 1}}), 
(b/.{x ^ x, y ^ y, z ^ z}/.{x ^{1,0, 0}, y ^ {0, 1, 0}, z ^ {0, 0, l}})].{x, y, z} 



Correcting a small glitch We will get some dodgy terms arising from the 
program attempting to also contract minus signs as separate entities. We can 
patch over this with the following rule: 

cross[_, — 1]:=0 

A. 3 Contracting polynomials now that the cross-products 
have been removed 

Now the cross products in contractions have been all re-expressed as simple 
sums of scalar products, we can continue with the contractions as normal. 

Taking summations outside The result from all the above simplifications 
might still be expandable: so any summations must therefore be taken outside 
the contraction sign again. 

C3[{al_.},Plus[x.,y._],{cl_.}]:= 

Module[{i}, xy = Expand[x + y]: Sum[C3[{al}, xy[[i]], {cl}], {i, 1, Length[xy]}]] 

Executing the contractions Now we need to find an efficient way to take 
each term and extract the xyz-depcndcnt components. 

First we separate out the terms containing xyz from the rest — this is so 
constant factors etc. don't get accidentally contracted when they shouldn't. 



Times[prefactor_/;FreeQ[FullForm[pref actor], x]/;FreeQ[FullForm[pref actor], y]/; 
FreeQ[FullForm[pref actor], z], xyz__], {cl__}]:= 



cl] 



Defining the cross product operator 



cross[a_, b_]: 



C3[{al_.} 
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Times[pref actor]C3[{al}, Times[xyz], {cl}]; 



Now we convert the polynomial (which is now in xyz form) into a tensorial 
form that can be used directly by the function Con. 

C3[{al__}, Times[xyzl_, xyz2__], {cl__}]:=Module[{i, power, n, xyz}, 
xyz:=List[xyzl,xyz2]/.x_P°"®''-- >Table[x, {i, 1, power}]; 
n = Length[xyz]; 

Con[T[al],T[Flatten[Table[xyz[[i]], {i. 1, n}]]], T[cl]] 
] 

C3[{al__}, Power[xyz_, power_], {cl_}]:=Module[{i, n}, 
Con[T[al],T[Table[xyz, {i, 1, power}]], T[cl]] 

] 

C3[al_,0,cl_]:=0; 

A. 4 Defining the scalar product symbofically 

Clear[scal, ^]; 

Attributes[scal] = {Orderless}; 

seal [vector., eij :=vector.{x, y, z}[[i]]; 

scal[vecl_, vec2_]:=vecl.vec2; 

A. 5 Application: chiral system 

The contraction mechanism has been defined to work on any kind of polynomial 
tensor involving a cross-product with another vector. Now we test it on the case 
of a chiral system, which has the field tensor 



Dipole-dipole correction We can see that this result is what you would ex- 
pect from a contraction of — with e' and e — the cross-product term in a 
has no effect as it has rank 3: 



Mixed (dipole-quadrupole) correction Here only the term in a affects 
the result. 




(44) 



C [{£'}, (x2-y2) + azzx(x2-y2),{e}] 



2e.xe'.x - 2e.ye'.y 



(45) 



C [{£'}, (x2-y2) + azzx(x2-y2),{k,e}] 



a (4(A:.ze.y + fc.ye.z)e'.x + 4(fc.ze.x + fc.xe.z)e'.y 
-|-4(A;.ye.x + /c.xe.y)e'.z) 



(46) 
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Quadrupole-quadrupole correction Now we are looking at tlie wliole ex- 
pression for tlic field tensor 

C [{k', e'}, (x2 - y2) + azz x (x^ - y^) , {k, e}] 

2fc.e'fc'.xe.x - 2fc.e'fc'.ye.y + 2A:.fc'e.xe'.x 
+2fc.x(/c'.xe.e' + fc'.ee'.x) - 2k.k'e.ye'.y - 2fc.y(fc'.ye.e' + fc'.ee'.y) ^ ' 
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